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ABSTRACT 

We perform numerical simulations of the Kelvin-Helmholtz instability in the mid-plane of a pro¬ 
toplanetary disk. A two-dimensional corotating slice in the azimuthal-vertical plane of the disk is 
considered where we include the Coriolis force and the radial advection of the Keplerian rotation flow. 
Dust grains, treated as individual particles, move under the influence of friction with the gas, while 
the gas is treated as a compressible fluid. The friction force from the dust grains on the gas leads to a 
vertical shear in the gas rotation velocity. As the particles settle around the mid-plane due to gravity, 
the shear increases, and eventually the flow becomes unstable to the Kelvin-Helmholtz instability. 
The Kelvin-Helmholtz turbulence saturates when the vertical settling of the dust is balanced by the 
turbulent diffusion away from the mid-plane. The azimuthally averaged state of the self-sustained 
Kelvin-Helmholtz turbulence is found to have a constant Richardson number in the region around the 
mid-plane where the dust-to-gas ratio is significant. Nevertheless the dust density has a strong non- 
axisymmetric component. We identify a powerful clumping mechanism, caused by the dependence of 
the rotation velocity of the dust grains on the dust-to-gas ratio, as the source of the non-axisymmetry. 
Our simulations confirm recent findings that the critical Richardson number for Kelvin-Helmholtz 
instability is around unity or larger, rather than the classical value of 1/4. 

Subject headings: diffusion — hydrodynamics — instabilities — planetary systems: protoplanetary 
disks — solar system: formation — turbulence 


1. INTRODUCTION 

One of the great unsolved problems of planet forma¬ 
tion is how to form planete simals, the kilo meter-sized 
precursors of real planets l|Safronovl 119691). At this 
size solid bodies in a protoplanetary disk can attract 
each other through gravitational two-body encounters, 
whereas gravity is insignificant between smaller bodies. 
Starting from micrometer-sized dust grains, the initial 
growth is c aused by the random Brownian motion of the 
grain s (e.gjBlum fc Wiirmll2?)nnl: [ Dullemond fc DominikI 
1200,^ see IHenning et al .1 1120061) for a review). The ver¬ 
tical component of the gravity from the central object 
causes the gas in the disk to be stratified with a higher 
pressure around the mid-plane. Even though the dust 
grains do not feel this pressure gradient, the strong fric¬ 
tional coupling with the gas prevents small grains from 
having any significant vertical motion relative to the gas. 
However, once the grains have coagulated to form peb¬ 
bles with sizes of a few centimeters, the solids are no 
longer completely coupled to the gas motion. They are 
thus free to fall, or sediment, towards the mid-plane of 
the disk. The increase in dust density opens a promising 
way of forming planetesimals by increasing the local dust 
density around the mid-plane of the disk to values high 
enough for gravitational fragmentation o f the dust layer 
jSafronovlll^ iGoldreich fc Wardl[T973|). 

There are however two major unresolved problems 
with the gravitational fragmentation scenario. Any 
global turbulence in the disk causes the dust grains 
to diffuse away from the mid-plane, and thus the 
dust density is kept at values that are too low for 
fragmentation. A turbulent a-value of 10“"^ is gen- 
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erally enough to pre vent efficient sedimentat i on to - 
wards the mid-plane l|Weidenschilling fc Cuzzil Il99,‘l|) . 
where as the a-value due to magnetorotational turbu¬ 
lence llB albus fc Hawlevlll991l: iBr andenburg et al.l[T99f)l: 
:mitagelll99 


iHawlev et al.lll995HArniitageili99fir is from a few times 

10“^ (found in local box simulations with no imposed 
magnetic field) to 0.1 and higher (in global disk simula¬ 
tions). The presenc e of a magnetica.llv dead zone around 
the disk mid-plane Gammigll996t iFromang et al.ll2?)?)3 


ISemenov et al.ll2nn4 ) may not mean that there is no tur¬ 

bulence in the mid-plane, as other instabili ties may set in 
and produce significant turb ulent motion llLi et al 11^ 
IKlahr fc Bodenheimeill^00,‘ll) . The magnetically active 
surface layers of the disk can even induce enough tur¬ 
bulent motion in the mid-plan e to possibly prevent ef - 
ficient sedimentation of dust l|Fleming fc Stonel 1200,‘ID . 
The presence of a dead zone may actually in itself be 
a source of turbulence. The sudden fall of the accre¬ 
tion rate can lead to a pile up of mass in the dead 
zone, poss ibly igniting the mag netorotational instability 
in bur sts llW iinsch et al.ll200fil)^ or a Rossby wave insta¬ 
bility llVarniere fc Taggeiil200,^ . 

The second major problem with the gravitational frag¬ 
mentation scenario is that even in the absence of global 
disk turbulence, the dust sedimentation may in itself 
destabilize the disk. Protoplanetary disks have a radial 
pressure gradient, because the temperature and the den¬ 
sity fall with increasing radial distance from the central 
object, so the gas rotates at a speed that is slightly below 
the Keplerian value. The dust grains feel only the gravity 
and want to rotate purely Keplerian. Glose to the equa¬ 
torial plane of the disk, where the sedimentation of dust 
has increased the dust-to-gas ratio to unity or higher, the 
gas is forced by the dust to orbit at a higher speed than 
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far away from the mid-plane where the rotation is still 
sub-Keplerian. Thus there is a vertical dependence of the 
gas rotation velocity. Such shear flow can be unstable 
to the Kelvin-Helmholtz instability (KHI), depending on 
the stabilizing effect of vertical gravity and density strat¬ 
ification. A necessary criterion for the KHI is that the 
energy required to lift a fluid parcel of gas and dust verti¬ 
cally upwards by an infinitesimal distance is available in 
the relat ive vertical motion be tween infinitesimally close 
parcels llChandrasekhail 1 196 111 . The turbulent motions 
resulting from the KHI are strong enough to puff up 
the dust layer and prevent the formation of an infinites- 
imally thin dust sheet around the mid-plane of the d isk 
llWeidenschillinelllQSntlAAidcnschilling fc CuzzilllQQ,*!! . 

Modifications to the gravitational fragmentation sce¬ 
nario have been suggested to overco me the problem of 
Kelvin-Helmholtz turbulence. iSekival (1998, hereafter re¬ 
ferred to as S98) found that if the mid-plane of the disk 
is in a state of constant Richardson number, as expected 
for small grains whose settling time is long compared 
to the growth rate of the KHI, then an increase in the 
global dust-to-gas ratio can lead to the formation of a 
high density dust cusp very close to the mid-plane of the 
disk, reaching potentially a dust-to-gas ratio of 100 al¬ 
ready at a global dust-to-gas ratio that is 10 times the 
canonical interstellar value of 0.01. The appearance of a 
superdense dust cusp in the very m id-plane has been in¬ 
terpreted bv lYoudin fc ShiJ ll2002ll as an inability of the 
gas (or of the KHI) to move more mass than its own away 
from the mid-plane. As a so urce of an increased va lue of 
the global dust-to-gas ratio, lYoudin fc Shul ll2002ll sug¬ 
gest that the dust grains falling radially inwards through 
the disk pile up in the inner disk. A slowly growing radial 
self-gravity mode in the dust density has also been sug¬ 
gested as the source of an increased dus t-to-gas ratio at 
certain radial locations l|Youdinll2005albr i. Trapping dust 
boulders in a turbulent flow is a mechanism for avoiding 
the proble m of self-induced Kelvin-He l mholtz turbulence 
altog ether tBaTge_^_Sonmieria| llM5l: jKjahrJ^Hennind 
19971 iHodgson fc Brandenburd Il998t iChavanisI I2QOOI: 
Johansen et al.ll2?)0^ . If the dust can undergo a gravi¬ 
tational fragmentation locally, because the boulders are 
trapped in features of the turbulent gas flow such as 
vortices or high-pressure regions, then there is no need 
for an extremely dense dust l a yer ar ound the mid-plane. 
iJohansen. Klahr. fc Hennind 11200(111 found that meter¬ 
sized dust boulders are temporarily trapped in regions of 
slight gas overdensity in magnetorotational turbulence, 
increasing the dust-to-gas ratio locally by up to two 
orders of magnitude. They estimate that the dust in 
such regions should have time to undergo gravitational 
fragme ntation before the h i gh-pre ssure regions dissolve 
again. iFromang fc Nelsor] J200,^ . on the other hand, 
find that vortices can even form in magnetorotationally 
turbulent disks, keeping dust boulders trapped for hun¬ 
dreds of disk rotation periods. The KHI cannot operate 
inside a vortex because there is no radial pressure gradi- 
ent, and thus no vertical shea r, in the center of the vortex 
llKlahr fc Bodenheimeill200(lll . 

From a numerical side it has been shown many times 
that a pure shear flow, i.e. one that is not explicitly sup¬ 
ported b^;_any_JorcesMsunstablej_both__wi^in^netic 
fields ((Keg^is^^I ]jl999UKeDDens fc T6thlll999ti and 
without l|Balbus et al.lll99m . But the key point here is 


that the vertical shear formed in a protoplanetary disk 
is due to the sedimentation of dust, and that the shear 
is able to regenerate as the dust falls down again, thus 
keeping the flow unstable to KHI. The description of the 
full non-linear outcome of such a system requires numer¬ 
ical simulations that include dust that can move relative 
to the gas. 

Linear stability analysis of dust-induced shear flows 
in protoplanetary disks hav e been performed fo r 
simplified physical conditions llSekiva fc IshitsiJ 1200011 . 
but also with Coriolis fo r ces a n d Ke plerian shear 
included lllshits u fc Seki^^ l2002l 1200,811 . Recently 
iGomez fc Ostrikeil ('2005. hereafter referred to as GO05) 
took an approach to include the dust into their numer¬ 
ical simulations of the Kelvin-Helmholtz instability by 
having the dust grains so extremely well-coupled to the 
gas that they always move with the instantaneous ve¬ 
locity of the gas. This is indeed a valid description of 
the dynamics of tiny dust grains. However, the strong 
coupling to the gas does not allow the dust grains to fall 
back towards the mid-plane. Thus the saturated state of 
the Kelvin-Helmholtz turbulence can not be reached this 
way. 

In this paper we present computer simulations where 
we have let the dust grains, represented by particles each 
with an individual velocity vector and position, move rel¬ 
ative to the gas. This allows us to obtain a state of self- 
sustained Kelvin-Helmholtz turbulence from which we 
can measure quantities such as the diffusion coefficient 
and the maximum dust density. A better knowledge of 
these important characteristics of Kelvin-Helmholtz tur¬ 
bulence is vital for our understanding of planet forma¬ 
tion. 

2. DYNAMICAL EQUATIONS 

We start by introducing the dynamical equations that 
we are going to solve for the gas and the particles. 

We consider a protoplanetary disk as a plane rotating 
with the Keplerian frequency l7o at a distance r = tq 
from a protostellar object. The plane is oriented so 
that only the azimuthal and vertical directions (which 
we name y and z, respectively) of the disk are treated. 
The absence of the radial a;-direction means that the Ke¬ 
plerian shear is ignored. The onset of the KHI is very 
likely to be affected by the presence of radial shear, since 
the unstable modes of the KHI are non-axisymmetric and 
therefore will be sheared out, but we believe the problem 
of Kelvin-Helmholtz turbulence in protoplanetary disks 
to be rich enough to allow for such a simplification as 
a first approach. There is of course the risk that the 
nature of the instability could c hange significantly wit h 
the inclusion of Keplerian shear ijlshitsii fc Sekiv'^l2nn.8ti . 
but on the other hand, the results that we present here 
regard mostly the dynamics of dust particles in Kelvin- 
Helmholtz turbulence, and we expect the qualitative re¬ 
sults to hold even with the inclusion of the Keplerian 
shear. 

As a dynamical solver we use the Pencil Code, a finite 
difference code that uses sixth order centered derivatives 
in space and a t hird order Runge-Ku tta time integra¬ 
tion scheme^. See iBrandenburel (1200.811 for details on the 
numerical schemes and test runs. 

^ The code, including modifications made for the current work, 
is available at 
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2 . 1 . Gas Equations 

The three components of the equation of motion of the 
gas are 

+ {u ■V)ux = 2noUy 

j (1) 

7 n 

duy , —s 1 „ 

+ (u ■ V)Uy = --S2oUj; 

1 dP e , , 

-^- [Uy-Wy], (2) 

P ay Tf 

duz , 1 dP e , X 

— + {u ■ V)uz = -- {Uz-Wz) . (3) 

at p dz Tf 

Here u = {ux,Uy,Uz) denotes the velocity field of the 
gas measured relative to the Keplerian velocity. We ex¬ 
plain now in some detail the terms that are present on 
the right-hand-side of equations 0 - 0 . The X- and y- 
components of the equation of motion contain the Cori¬ 
olis force due to the rotating disk, to ensure that there 
is angular momentum conservation. Since velocities are 
measured relative to the Keplerian shear flow, the advec- 
tion of the rotation flow by the radial velocity component 
has been added to the azimuthal component of the Cori¬ 
olis force, changing the factor —2 to —1/2 in equation 
0 . We consider local pressure gradient forces only in 
the azimuthal and vertical directions, whereas there is a 
constant global pressure gradient force in the radial di¬ 
rection. The global density is assumed to fall radially 
outwards as din p/9 In r = a, where a is a constant. As¬ 
suming for simplicity that the density decreases isother- 
mally, we can write the radial pressure gradient force as 


1 dP 
p dr 


1 2 In p 

7 


dr 


( 4 ) 


Here 7 = 5/3 is the ratio of specihc heats and Cs is the 
constant sound speed. Rewriting this expression and us¬ 
ing the isothermal disk expression H = Cs/i7oj we arrive 
at the expression 


p dr 


1 iJ 9 In p 




( 5 ) 


7 r 91nr 

We then proceed by defining the dimensionless disk pa¬ 
rameter 9 = -f- , where H/r is the scale-height to ra¬ 

dius ratio of the disk. This parameter can be assumed to 
be a constant for a protoplanetary disk. Using equation 
© and the definition of (3 leads to the global pressure 
gradient term in equation ©. We use throughout this 
work a value of /3 = —0.1. 

The ratio between the pressure gradient force Ap and 
two times the solar gravity, 


V = 


^9 _ p dr 

2 g 


1 dP 

p dr 
212 gr 


( 6 ) 


is often used to pa rameterize sub-Keplerian disks 
l|Nakagawa et Assuming again an isothermally 

falling density, equation can be written as 


1 1 ( H 

^ = -2-A~ 


91np 
91nr ’ 


( 7 ) 


http://www.nordita.dk/software/pencil-code/ 


The connection between our (3 and the more widely used 
77 is then 


V = - 7 ; -/?• 

z 7 r 


( 8 ) 


The last term in equations 0-0 is the fric¬ 
tion force that the dust particles exert on the gas. 
We discuss this in further detail in Sect. 12.31 be¬ 
low. Here the dust velocity field w — {wx,Wy,Wz) 
is a map of the particle velocities onto the grid. 
To stabilize the finite difference numerical scheme of 
the Pencil Code, we use a s ixth-order momentum- 
conserving hyperviscosity (e.g. |Bra ndeiibuig_£^_Sarsoi] 
2^22 ; iHauaen fc Brandenburel 12004 l.lohansen fc Klahr 
200^. Hyperviscosity has the advantage over regular 
second-order viscosity in that it has a huge effect on un¬ 
stable modes at the smallest scales of the simulation, 
while leaving the largest scales virtually untouched. 


2 . 2 . Mass and Energy Conservation 

The conservation of mass, given by the logarithmic 
density Inp, and entropy, s, is expressed in the conti¬ 
nuity equation and the heat equation, 

^ + (m-V) lnp = -V-w, (9) 

ot 

|^ + (m-V)s = 0. (10) 

The advection of the global density gradient due to any 
radial velocity has been ignored, as well as viscous heat¬ 
ing of the gas. We calculate the pressure gradient force 
in equations 0-0 by rewriting the vector term as 


-p-iVP = -c2(Vs/cp + Vlnp) (11) 

and using the ideal gas law expression 


2^2 
Cs = 7 — = Cso exp 
P 


7 s/cp -I- ( 7 -I) In ■ 


PoJ 


( 12 ) 


The two constants Cso and po are integration constants 
from the integration of the first law of thermodynam¬ 
ics. We have chosen the integration constants such that 
s = 0 when Cg = Cso and p = po. To allow for grav¬ 
ity waves, we must use the perfect gas law, rather than 
a simple polytropic equation of state. We stabilize the 
continuity equation and the entr opy equation by usi ng a 
5th order upwinding scheme fsee lDobler et alJl 200 (TH for 
the advection terms in equations 0 and m- 


2.3. Dust Equations 

The dust grains are treated as individual particles mov¬ 
ing on the top of the grid. Therefore they have no advec¬ 
tion term in their equation of motion, whose components 
are 



The index i runs in the interval i = 1... A, where N 
is the number of particles that are considered. The last 
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TABLE 1 
Run parameters 


Run 

l?0Tf 

eo 

Ly 

X 

L. 

Ny 

X 

N, 

N 

Ni 


A 

0.02 

0.01 

0.40 

X 

0.20 

256 

X 

128 

400,000 

12.2 

200 

B 

0.10 

0.01 

0.40 

X 

0.20 

256 

X 

128 

400,000 

12.2 

200 

C 

1.00 

0.01 

0.10 

X 

0.05 

256 

X 

128 

400,000 

12.2 

80 

Be2 

0.10 

0.02 

0.40 

X 

0.20 

256 

X 

128 

400,000 

12.2 

100 

Be5 

0.10 

0.05 

0.40 

X 

0.20 

256 

X 

128 

400,000 

12.2 

100 

BelO 

0.10 

0.10 

0.40 

X 

0.20 

256 

X 

128 

400,000 

12.2 

100 

Br512 

0.10 

0.01 

0.40 

X 

0.20 

512 

X 

256 

1,600,000 

12.2 

100 


Note. — Col. (1): Name of run. Col. (2): Stokes number. Col. (3): 
Global dust-to-gas ratio. Col. (4): Size of simulation box. Col. (5): Grid 
resolution. Col. (6): Number of particles. Col. (7): Number of particles per 
grid point. Col. (8): Total time of run in units of 


terms in equations is the friction force. The 

friction force is assumed to be proportional to the veloc¬ 
ity difference between dust and gas with a characteristic 
braking-down time-scale of Tf, called the friction time. 
To conserve the total momentum, the dust must affect 
the gas by an oppositely directed friction force with fric¬ 
tion time Tf/e, as included in the last terms of equa¬ 
tions Q-©- Here e is the local dust-to-gas mass ra¬ 
tio Pd/p- The dust density pd at a grid point is cal¬ 
culated by counting the number of particles within a 
grid cell volume around the point and multiplying by 
the mass density pd that each particle represents. The 
mass density per particle depends on the number of par¬ 
ticles and on the assumed average dust-to-gas ratio cq as 
Pd = eop/-^i, where is the number of particles per 
grid cell. Since the gas is approximately isodense and 
isothermal, we can assume that the friction time is inde¬ 
pendent of the local state of the gas at the position of 
a particle. We also assume that the friction time does 
not depend on the velocity difference between the par¬ 
ticle and the gas. This is valid in the Epstein regime, 
but also in the St okes regime when the flow around the 
grains is laminar llWeidenschillinglll977ll . For conditions 
typical for a protoplanetary disk at a radial distance of 5 
AU from the central object, a given Stokes number 
corresponds to the gr ain radius measured in meters (e.g. 
iJohansen et al.ir 200 (lll . although this depends somewhat 
on the adopted disk model. We include the friction force 
contribution to the computational time-step 5t by requir¬ 
ing that the time (Jtjfric = Tf /(1 -|- e) is resolved at least 
five times in a time-step. This restriction is strongest for 
small grains and for large dust-to-gas ratios, whereas the 
Courant time-step of the gas dominates otherwise. 

The particle positions change due to the velocity of the 
particles as 


dt 

dt 

dt 


= 0 , 

(16) 

= 4 ^ 

(17) 

= u«. 

(18) 


Because the simulations are done in two dimensions, we 
have not allowed particles to move in the cc-direction. 
The particles are still allowed to have a radial velocity 
component. This is equivalent to assuming that all ra¬ 


dial derivatives are zero, so that no advective transport 
occurs in this direction. 


3. RICHARDSON NUMBER 


Before we discuss the setup of the numerical simula¬ 
tions and the results, we describe in this section some of 
the analytical results that are already known about the 
KHI. 

The stability of a shear flow can be characterized 
through the Richardson number Ri, defined as 


Ri = 


gzd\n{p + pd)/dz 

{duyjdzY 


(19) 


The Richardson number quantifies the fact that verti¬ 
cal gravity and density stratification of gas and dust 
d\n{p -\- pd)/dz are stabilizing effects, whereas the shear 
duyjdz is destabilizing. As shown by Chandrasekhar 
from very simple considerations of the free energy that is 
present in a stratified shear flow, a flow with Ri > 1/4 is 
always stable, whereas Ri < 1/4 is necessary, bu t not suf¬ 
ficient, for an instability l|Chandrasekhailfl96lL p. 491). 
These derivations do, however, not include the effect of 
the Coriolis force, a point that we shall return to later. 

For dust-induced shear flows in protoplanetary disks, 
S98 derived an expression for the vertical density distri¬ 
bution of dust in a protoplanetary disk that is marginally 
stable to the KHI, i.e. where the gas flow has a constant 
Richardson number equal to the critical Richardson num¬ 
ber for stability Ric. For small dust grains the dust-to- 
gas ratio e(z) in this state can be written as 


e{z) 


0 


1 for \z\ < Zd 
for \z\ > Zd 


( 20 ) 


where ei is the dust-to-gas ratio in the mid-plane, Zd = 
Hd^l-1/{1 -l-ei)^, and Hd is the width of the dust 
layer. The effect of self-gravity between the dust grains 
has been ignored. The width of the dust layer can fur¬ 
thermore be written as 


Hd = V^J-^H, ( 21 ) 

27 

where (3 is the radial pressure gradient parameter intro¬ 
duced in Sect. im 7 is the ratio of specific heats and H 
is the scale-height of the gas. For Ric = 1/4 and 7 = 5/3 

Hd/H=^\(3\, (22) 



















Dust Sedimentation and Kelvin-Helmholtz Turbulence 


5 



Fig. 1.— The onset of the Kelvin-Helmholtz instability for cm-sized pebbles with = 0.02. The initial Gaussian particle distribution 

falls towards the mid-plane of the disk on the characteristic time-scale of fgrav = ^ . The increased vertical shear in the gas 

rotation velocity eventually makes the disk unstable to the KHI, forming waves that finally break as the turbulence goes into its non-linear 
state. 


SO the width of the marginally stable dust layer is a few 
percent of a gas scale height. 

The expression in equation GDI) allows for two types 
of dust stratihcation in the marginally stable disk. For 
Pd ^ P the dust density is constant around the mid¬ 
plane, whereas for pd ^ p, a cusp of very high dust 
density can exist around the mid-plane (S98). Such a 
cusp can form for two reasons when the dust-to-gas ratio 
is above unity. Firstly because the gas flow is forced 
to be Keplerian in such a large region around the mid¬ 
plane that the vertical shear is reduced there, stabilizing 
against the KHI, and secondly because it requires a lot of 
energy to lift up so much dust awa y from the mid - plane . 
This effect has been interpreted bv lYoudin fc ShiJ 11200211 
as the gas only being able to lift up its own equivalent 
mass due to KHI. 

4. INITIAL CONDITION 


The initial condition of the gas is an isothermal, strat¬ 
ified disk with a scale height H. The density depends on 
the height over the mid-plane ^ as 

p(.) = ^ (23) 

where pi is the density in the mid-plane. The scale height 
is H = Cs/I2oj where Cg determines the constant initial 
temperature, to sustain hydrostatic equilibrium in the 
vertical direction. From the definition of the gas col¬ 
umn density S, we can calculate the mid-plane density 
as Pi = There is no similar equilibrium to 

set the dust density in the disk. Thus we distribute the 
particles in a Gaussian way around the mid-plane with 
a scale height Hd, a free parameter, and normalize the 
distribution so that Ed = ^qE, where Cq is the global 
dust-to-gas ratio in the disk. 

The constant global pressure gradient force effectively 
decreases the radial gravity felt by the gas, and thus the 
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orbital speed is no longer Keplerian, but slightly smaller. 
The sub-Keplerian velocity is given by 

„( 0 ) = A 


27 


(24) 


This expression is obtained by setting Ux = dux/dt = 
e = 0 in equation CJ. The dust grains, on the other 
hand, do not feel the global pressure gradient and would 
thus in the absence of friction move on a Keplerian orbit 
with Vy'^ — 0. The drag force from the gas, however, 
forces the dust grains to move at a speed that is below 
the Keplerian value, at least when the dust-to-gas ratio 
is low. When the dust-to-gas ratio approaches unity or 
even larger, the gas is forced by the dust to move with 
Keplerian speed. The equilibrium gas and dust veloc- 
ity can be calculated a s a function of dust-to-gas ratio 
l|Na,ka,gawa, et but we choose to simply start 

the gas with a sub-Keplerian velocity and the dust with 
a Keplerian velocity, and then let the velocities approach 
the equilibrium dynamically (this happens within a few 
friction times). This way we have checked that the nu- 
merical solution indee d approaches the expressions by 
iNa.ka.gawa, et a 1.1 ill 98fifl for all the velocity components 
of the gas and the dust, which serves as a control that 
the friction force from the dust particles on the gas has 
been correctly implemented in the code. 

In the equilibrium state the gas has a positive radial 
velocity in the mid-plane, but this does not lead to any 
change of the gas density, since we have ignored the ad- 
vection of the global density in the continuity equation. 
The effect of an outwards-moving mid-plane on the global 
dynamics of a protoplanetary disk is a promising topic of 
future research, but it is beyond the scope of this paper. 

Periodic boundary conditions are used for all variables 
in the azimuthal direction. In the vertical direction we 
impose a condition of zero vertical velocity, whereas the 
two other velocity components have a zero derivative con¬ 
dition over the boundary. The logarithmic mass den¬ 
sity and the entropy have a condition of vanishing third 
derivatives over the vertical boundary. 

We run simulations for three different grain sizes, 
HoTf = 0.02,0.1,1.0, respectively. When assuming com¬ 
pact spherical grains, these numbers correspond to sizes 
of centimeters (pebbles), decimeters (rocks), and meters 
(boulders), respectively. The computation parameters 
are listed in Tabled The size of the box is set accord¬ 
ing to the vertical extent of the dust layer in the state 
of self-sustained Kelvin-Helmholtz turbulence. We make 
sure that the full width of the dust layer fits at least five 
times vertically in the box to avoid any effect of the ver¬ 
tical boundaries on the mid-plane. The azimuthal extent 
of the box is set so that the full width of the dust layer 
fits at least ten times in this direction. Thus the unstable 
modes of the Kelvin-Helmholtz turbulence, which have a 
wavelength that is similar to the width of the dust layer, 
are very well-resolved. 

Since the ratio of particles to grid points is « 12, 
the computation time is strongly dominated by the par¬ 
ticles. Each particle needs to “know” the gas velocity 
field at its own position in space, to calculate the fric¬ 
tion forces. For parallel runs we distribute the particles 
among the different processors so that the position of 
each particle is within its “host” processor’s space in¬ 
terval. As we show in Sect. |3 the particles tend to 
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Fig. 2. — Contour plot of the dust density of cm-sized pebbles 
averaged over the azimuthal y-direction, as a function of time t and 
height over the mid-plane 2 . After the Kelvin-Helmholtz instabil¬ 
ity sets in and saturates into turbulence, the width of the dust 
layer stays approximately constant. The black regions contain no 
particles at all. 
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Fig. 3. — Same as Fig.|21 but for dm-sized rocks with SloTf = 0.1. 
The sedimentation time-scale is much faster than in Fig.|^ but the 
width of the dust layer in the self-sustained state of turbulence is 
approximately the same. 



have a strongly non-axisymmetric density distribution in 
the Kelvin-Helmholtz turbulence. This clumping means 
that the number of particles on the individual proces¬ 
sors varies by a factor of around five, thus slowing the 
code down by a similar factor compared to a run where 
the particles were equally distributed over the processors. 
For this reason we used 32 Opteron processors, each with 
2.2 GHz CPU speed and Infiniband interconnections, for 
around one week for each run. 
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Fig. 4.— Dust density and Richardson number of grains with 
= 0.02 averaged over the azimuthal direction and over time. 
The dust-to-gas ratio in the mid-plane is close to unity and falls 
rapidly outwards. The Richardson number is approximately con¬ 
stant in the mid-plane and has a value around unity. 
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zlH 

Fig. 5.— Same as in Fig. but for HoTf = 0.1. Although the 
density in the mid-plane is similar to the value for smaller grains, 
the dust is more settled and has less pronounced wings away from 
the mid-plane. The Richardson number is again around unity in 
the mid-plane. 


5. DYNAMICS AND DENSITY OF SOLIDS 

In this section we focus on the dynamics and the den¬ 
sity of the dust particles. The linear growth rate of the 
Kelvin-Helmholtz instability and the statistical proper¬ 
ties of the Kelvin-Helmholtz turbulence are treated in 
the next two sections. 

Some representative snapshots of the particle positions 
for run A (I^oFf = 0.02, or cm-sized pebbles) are shown 
in Fig. ^ The particles, with an initial Gaussian den¬ 
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Fig. 6.— Same as Fig.but for m-sized boulders with f2oR = 
1.0. The equilibrium scale height of the dust is around 10 times 
lower than for the smaller grains. 



ZlH 



ZlH 

Fig. 7.— Same as in Fig.|^ but for Corf = 1.0. The Richardson 
number around the mid-plane is here way lower than for the smaller 
grains, around 0.1. This is caused by the extremely rapid settling 
of m-sized boulders to the mid-plane. 


sity distribution, settle to the mid-plane due to gravity, 
on the characteristic time-scale tgrav = l/('R^o) ~ ^0- 
When the width of the layer has decreased to around 
0.01 scale heights (the two middle panels of Fig. ^1, some 
wave pattern can already be seen in the dust density. It 
is the most unstable Uz{y) mode, with a wavelength that 
is comparable to the vertical width of the layer, that is 
growing in amplitude. Some shear times later, in the 
two bottom panels of Fig. ^ the growing mode enters 
the non-linear regime, and impressive patterns of break¬ 
ing waves appear. The simulation then goes into a state 
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Fig. 8.— Contour plot of the particle density for Corf = 1.0. A clump, indicated by the arrow, oscillates around the mid-plane, a type 
of motion that is only allowed for particles with CpTf > 0.5. 


of fully developed Kelvin-Helmholtz turbulence. 

5.1. Pebbles and Rocks 

In Figs. 12 and |2 we show azimuthally averaged dust 
density contours as a function of time for grains with 
friction time IJoTf = 0.02 (run A) and UoTf = 0.1 (run 
B, decimeter-sized rocks), respectively. These two grain 
sizes show very similar behavior with time. In the be¬ 
ginning the particles move towards the mid-plane unhin¬ 
dered because of the lack of turbulence. The sedimenta¬ 
tion happens much faster in run B than in run A because 
of the different grain sizes. When the KHI eventually sets 
in, the dust layer is puffed up again and quickly reaches 
an equilibrium configuration where the vertical distribu¬ 
tion of dust density is practically unchanged with time. 

One can already here suspect that the equilibrium dust 
density is indeed, as predicted analytically by S98, dis¬ 
tributed in such a way that the flow has a constant 
Richardson number. In Figs. 0] and El we plot the time- 
averaged dust density and the Richardson number as a 
function of vertical height over the mid-plane, again for 
the two small grain sizes. The dust-to-gas ratio reaches 
unity in the mid-plane and drops down rapidly away 
from the mid-plane. For run A the Richardson num¬ 
ber is approximately constant, just above unity, in the 
region around the mid-plane that has a significant dust 
density. For run B the value of the Richardson number 
is also constant, although somewhat smaller than for the 
centimeter-sized grains, because the more rapid sedimen¬ 
tation of these larger grains allows the disk to sustain a 
stronger vertical shear. 

5.2. Boulders 


For bodies with = 1-0 (run C, m-sized boulders), 
the azimuthally averaged dust density is shown in Fig. El 
These meter-sized boulders fall rapidly to the mid-plane, 
on a time-scale of one shear time, because they are not 
as coupled to the gas as smaller grains. The scale height 
of the boulders is very small, less than one percent of the 
scale height of the gas, because the grains are falling so 
fast that the disk can sustain a much lower Richardson 
number than the critical. This is also evident from Fig. [3 
The Richardson number is well below unity, around 0.1, 
where significant amounts of dust is present. 

A major difference between the large grains and the 
small grains is the presence of bands in Fig. El The dust 
grains are no longer smoothly distributed over z, but 
rather appear as clumps that oscillate around the mid¬ 
plane. The oscillation of a single clump is evident from 
Fig.El Here the dust density contours are plotted at four 
times separated by one shear time. The clump indicated 
by the arrow is oscillating around the mid-plane. Such 
oscillatory behavior is also be expected from the follow¬ 
ing considerations. Friction ensures that small grains ar¬ 
rive at the mid-plane with zero residual velocity, whereas 
larger grains perform damped oscillations around z = 0 
with a damping time of one friction time. The distinction 
between the two size regimes can be derived by looking 
at the differential equation governing vertical settling of 
particles. 


dvzjt) 

dt 


= —I2qZ- 1 

Ti 


(25) 


This second order, linear o rdinary differential eq uation 
can be solved trivially le.g. iNakagawa et alJIlQSfil) . The 
result is a split between two types of solutions, depend- 
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Fig. 9.— The root-mean-square 2 -coordinate of the particles for 
the different runs. For = 0.02 and = 0.1, the width 

of the dust layer is around 1% of a gas scale height, whereas for 
meter-sized boulders with Qqti = 1.0, the strong sedimentation 
results in a much lower width, only around 0.25% of a scale height. 


ing on the value of For l7o'T < 0.5 the solution is 

a purely exponentially decaying function. On the other 
hand, for > 0.5 the solutions are damped oscilla¬ 
tions with a characteristic damping time of around one 
friction time. For a friction time around unity, the am¬ 
plitude of the dust density in a laminar disk would still 
become virtually zero in just a few friction times. This 
is not the case in Fig. El where the dust scale height stays 
approximately constant for at least 80 friction times (the 
end of the simulation). The KHI is continuously pump¬ 
ing energy into the dust layer at the same rate as the 
oscillations are damped. 

In Fig. M we plot the root-mean-square value of the 
z-coordinate of the particles as a function of time for 
all three runs. The calculated “scale height” of the cm- 
sized and dm-sized particles are very similar, although 
the larger particles have a bit lower scale height than 
the smaller particles. The m-sized particles have a scale 
height of about one quarter of a percent of that of the 
gas. 


5.3. Maximum Density and Clumping 

It is of great relevance for planetesimal formation to 
find the highest dust density that is permitted in the 
saturated Kelvin-Helmholz turbulence. Both coagulation 
and gravitational fragmentation depend strongly on the 
mass density of the dust layer. High densities can occur 
not only when the gas flow or the size of the boulders 
allow for a high mid-plane density, but also in certain 
points of the turbulent flow where dust particles tend to 
accumulate. The latter can only be explored in computer 
simulations, so we examine the maximum dust density in 
any grid point in more detail in this section. 

The maximum mass density of dust particles in any 
grid cell is plotted in Fig. as a function of time. Even 
though the mid-plane dust-to-gas ratio is on the aver¬ 



tlQ.Q 

Fig. 10.— The maximum dust density in any grid cell as a 
function of time, in units of the gas mid-plane density, as a function 
of time. The value is much higher than the azimuthally averaged 
mid-plane densities (presented in previous figures). 


age of the order unity for all grain sizes, the maximum 
density is much higher at all times, especially for meter¬ 
sized particles where the maximum dust-to-gas ratio can 
be up to one thousand. Decimeter-sized particles have 
a maximum dust-to-gas ratio of around 20 at all times, 
whereas the value for centimeter-sized particles is around 
10. This is potentially important for building planetesi- 
mals. Even if the critical density for gravity-aided plan¬ 
etesimal formation is not reached globally, this is still 
possible in certain regions of the turbulent flow. Such a 
gr avoturbulent formation of planetesimals was proposed 
by l.lohansen et al.1 1120061) to lead to the formation of 
planetesimals in a magnetorotationally turbulent gas. 

In Fig. CH we plot contours of the vertically averaged 
dust density as a function of azimuthal coordinate y and 
time t for decimeter-sized rocks (run B). It is evident that 
the dust density has a strong non-axisymmetric compo¬ 
nent once the Kelvin-Helmholtz turbulence is fully de¬ 
veloped. Dense clumps are seen as white stripes, while 
regions of lower dust-to-gas ratio are gray. A simple way 
to quantify the amount of non-axisymmetry is to look at 
the mean deviation of the dust density from the average 
density. We define the azimuthal clumping factor Cy as 


Viinyiy) - {ny{yW) 

My)) 


(26) 


Here ny{y) = {n(y,z))z is the dust number density av¬ 
eraged over the vertical direction. Axisymmetry implies 
Cy = 0, whereas higher values of Cy imply stronger and 
stronger non-axisymmetry. We plot in Fig. 1121 the az¬ 
imuthal clumping factor as a function of time for all three 
values of the friction time. For centimeter and decimeter 
particles the clumping is of the order unity, or in other 
words, the average grid point has a density variation from 
the average that is on the same order as the average, i.e. 
very strong clumping. For meter-sized particles the az¬ 
imuthal clumping is even stronger. 

The tendency for the dust particles to clump is a con- 
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Fig. 11.— Vertically averaged dust density of rocks with 
l^oTf = 0.1 as a function of azimuthal coordinate y and time t. 
The clumping mechanism is evident from the plot. Regions of high 
dust-to-gas ratio (light) move slower than regions of low dust-to¬ 
gas ratio (gray), seen in the different slopes of the bright and dark 
wisps on the plot, causing the high density clumps to be contin¬ 
uously fed by low density material. One also sees the rarefaction 
tail going to the left of the dense clumps and the shock front that 
is formed against the sub-Keplerian stream. 


sequence of the sub-Keplerian velocity of the gas. Turn¬ 
ing again to Fig. CH one sees that brighter regions move 
at a lower speed (relative to the Keplerian speed) than 
dark regions do. The speed of a clump is evident from 
the absolute value of the angle between the tilted time- 
space wisp and the time-axis. Bright wisps have a higher 
angle with the time-axis than dark wisps. This insta- 
bil ity is very related to the s treaming instability found 
bvlYoudin fc Goodmaiil ll2005ll. although in our simula¬ 
tions the clumping happens in the {y, 2 ;)-plane and not 
the (a;,z)-plane as in the analysis by Youdin & Good¬ 
man. Still the instability is powered in both planes by 
the dependence of the dust velocity on the dust-to-gas 
ratio, so we consider the instability in the {y, z)-plane 



Fig. 12.— Azimuthal clumping factor Cy versus time for all 
three grain sizes. A value of around unity corresponds to strong 
clumping with the average point being 100% away in density from 
the average density. 


a special case of the stream ing instability. We refer to 
lYoudin fc Goodmanl ll200,^ for a linear stability analy¬ 
sis of the coupled gas and dust flow. In the rest of this 
section we instead focus on describing the non-linear out¬ 
come of the streaming instability with a simple analogy 
to a hydrodynamical shock. 

Using the derivations given bv iNakagawa et all lll98fill 
for the equilibrium gas and dust velocities as a function 
of the dust-to-gas mass ratio e, one can write up the 
azimuthal velocity component of the dust as 


Wy 


l + e 

(1 -I- e)2 -I- 


,(o) 


(27) 


Thus clumps with a high dust-to-gas ratio move slower, 
relative to the Keplerian speed, than clumps with a low 
dust-to-gas ratio. The small clumps crash into the big 
clumps and form larger structures. At the same time, the 
large clumps steepen up against the direction of the sub- 
Keplerian flow and develop an escaping tail downstream. 
This is qualitatively similar to a shock. Considering the 
continuity equation of the dust-to-gas ratio 


de de dwy 

dt dy ^ dy 


(28) 


and inserting equation m in the limit of small Stokes 
numbers l7orf ^ 1, equation can be reduced to 


de _ de 

dt (1 -I- e)2 dy ’ 


(29) 


qualitatively similar to the advection equation of fluid 
dynamics. The shock behavior of the clumps arises be¬ 
cause the advection velocity /{I ^-e)^ depends on the 
dust-to-gas ratio. 


5.4. Varying the Global Dust-to-gas Ratio 
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It is of great interest to investigate the dependence 
of the mid-plane dust density on the global dust-to-gas 
ratio in the saturated state of Kelvin-Helmholtz turbu¬ 
lence. Increasing the dust-to-gas ratio beyond the inter¬ 
stellar value should potentially lead to the creation of a 
very dense mid-plane of dust that the gas is not able to 
lift up, making the dust lay er dense enough to undergo 
gravitational fragmentation ll^Sekivalll998t lYoiidin fc Shd 
l2?inl lYoudin ChiandEifMir 

The analytically predicted mid-plane dust-to-gas ratio 
ei is found by applying the normalization 


pit{z)Az = Yd 


(30) 


to the constant Richardson number dust density of equa¬ 
tion (HD). Here Yd is the dust column density, a free pa¬ 
rameter that we set through the global dust-to-gas ratio 
eo as Yd = eo^- We have approximated the gas den¬ 
sity by the gas density in the mid-plane pi, because for 
z H, the variation in gas density is insignificant com¬ 
pared to the variation in dust density. We proceed by 
inserting the expression for the dust density in a disk 
with a constant Richardson number, from equation (EDI), 
into the integral in equation Defining the parame¬ 
ter _ 


X = 


y/ei{2 + ei) 

1 + ei 


(31) 


the integration yields 


- 2x -I- In 


1 + X 
1 -X 




(32) 


This is a transcendental equation that we solve numer¬ 
ically for X as a function of the input parameters i/d) 
given by equation ED, and Yd. Once x is calculated, 
then the dust-to-gas ratio in the mid-plane ei is known 
from equation ED. 

In Fig. El we plot the analytical mid-plane dust-to-gas 
ratio Cl as a function of the global dust-to-gas ratio eo 
(dotted line). The non-linear behavior of ei is evident, 
and already for eo = 0.1 does the mid-plane dust-to-gas 
ratio approach ei = 100, which should be enough to have 
a gravitational instability in the dust layer. We also run 
numerical simulations with an increased global dust-to- 
gas ratio (runs Be2, Be5 and BelO, see Table^l to see if a 
mid-plane dust density cusp develops as predicted. The 
resulting mid-plane dust-to-gas ratio is indicated with 
stars in Fig. El To avoid having a very low time-step, 
because of the strong friction that the dust exerts on the 
gas when the global dust-to-gas ratio is increased, we 
have locally increased the friction time in regions of very 
high dust-to-gas ratio. This approach conserves total 
momentum because the friction force on the gas and on 
the dust are made lower at the same time. In the regions 
where the friction time is increased, the dust-to-gas ratio 
is so high that gas is dragged along with the particles 
anyway, so the precise value of the friction time does not 
matter. As seen in Fig. El the mid-plane dust-to-gas 
ratio does indeed increase non-linearly with global dust- 
to-gas ratio, following a curve that is within a factor of 
two of the analytical curve. This gives support to the 
theory that an increased global dust-to-gas ratio, e.g. 
due to solids that are transported from the outer part 
of the disk into the inner part, can lead to such a high 
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Fig. 13.— Mid-plane dust-to-gas ratio ei as a function of global 
dust-to-gas ratio eo- The dotted line shows the analytical value for 
a dust density with a constant Richardson number of Ric = 1.0, 
while the stars show the result of the numerical simulations for 
different values of the global dust-to-gas ratio. The results agree 
nicely, giving support to the idea that a dust-to-gas ratio that is 
higher than the interstellar value can give rise to high enough mid¬ 
plane dust density for a gravitational instability in the dust layer. 
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dust-to-gas ratio in the disk mid-p lane that the dust lay er 
fragments to form planetesimals lIYoudin fc Shull20iilh . 


6. GROWTH RATES 


The simulations presented so far all imply a critical 
Richardson number that is of the order unity, rather 
than the classically adopted value of 1/4. This confirms 
the findings of GO05 that shear flows with a Richardson 
number that is higher than the classical value are actu¬ 
ally unstable when the Coriolis force is included in the 
calculations. To quantify the linear growth of the insta¬ 
bility we have run simulations of initial conditions with 
a constant Richardson number and measured the growth 
rates of the KHI. Because we are only interested in the 
linear regime, we have chosen for simplicity to treat dust 
as a fluid rather than as particles. Thus we solve equa¬ 
tions similar to equations for the dust velocity 

field w including an advection term {w ■ V)ia. A conti¬ 
nuity equation similar to equation m for the logarithmic 
dust number density Inn is solved at the same time. 

We consider initial conditions with a constant Richard¬ 
son number, in the range between 0.1 and 2.0, as given 
by equation ED- The dust-to-gas ratio in the mid-plane 
ei is known from equations ED and ED- To avoid 
any effects of dust settling, we set the friction time to 
iloTf = 0.001. The gravitational settling time is then as 
high as 1000f7(/^, and since this is much longer than the 
duration of the linear growth, the effect on the measured 
growth rate is insignificant. The vertical velocity of the 
dust is set to the terminal settling velocity Wz = —TfilgZ. 
Since the dust velocity is not zero, the gas will feel some 
friction from the falling dust. The vertical component of 
the equation of motion of the gas is 


duz 


-12/ 


1 2 ^ 111 P 


— {Uz - Wz) . 
Tf 


(33) 


We insert the terminal velocity expression for the dust 
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velocity into equation and look for equilibrium so¬ 
lutions with Uz = duz/dt = 0. The equation is then 
reduced to 

Q = -{l + e)nlz-^cl^. (34) 

The drag force exerted by the falling dust on the gas 
mimics a vertical gravity, and therefore we have com¬ 
bined it with the regular gravity term. In a way it is 
the gravity on the dust that the gas feels, only it is 
transferred to the gas component through the drag force. 
One can int erpret this as the gas feeling a stronger grav¬ 
ity I7g = \/l + el7o in places of high dust-to-gas ratio, 
which leads to the creation of a small cusp in the gas 
density around the mid-plane. Inserting now equation 
CT into equation m and applying the boundary con¬ 
dition p{z = 0) = Pi yields 

1 _ / I 1 

i+<=i IT+ttf. 

for \z\ < Zd 

_1 

2 IF FT+FF. 

for \z\ > Zd 

(35) 

The cusp around the mid-plane, caused by the extra 
gravity imposed by the falling dust on the gas, is shown in 
Fig. 1141 The variation in density from a normal isother¬ 
mal disk is only a few parts in ten thousand, so the effect 
is not big. On the other hand, it is important to have a 
complete equilibrium solution as the initial condition for 
the measurement of the linear growth of the instability, 
as otherwise dynamical effects could dominate over the 
growth. 

We measure the linear growth rate of the KHI by pre¬ 
scribing a dust-to-gas ratio according to equation TO 
and a gas density according to equation We then 

set the velocity field s of gas and dust a ccord ing to the ex¬ 
pressions derived in INakagawa et gdl lll98(ill . The width 
of a dust layer with a constant Ri depends on the value 
of Ri according to equation so we have made sure to 
always resolve the unstable wavelengths by making the 
box larger with increasing Ri. The fluid simulations are 
all done with a grid resolution of Ny x Nz = 256 x 128. 

The measured growth rates are shown in Fig. 1151 At a 
Richardson number close to zero, the growth rate is simi¬ 
lar in magnitude to the rotation frequency Uq of the disk, 
whereas for larger values of the Richardson number, the 
growth rate falls rapidly. We find that there is growth 
out to at least Ri = 2.0, with a growth rate approaching 
CO = O.OII^Q. There is no evidence for a cut-off in the 
growth rate at the classical value of the critical Richard¬ 
son number of 0.25. The range of unstable Richard¬ 
son numbers is in good agreement with the mid-plane 
Richardson number in the particle simulations shown in 
Figs. EJandEl This is another confirmation that the crit¬ 
ical Richardson number is around unity or higher when 
the Coriolis force is included in the calculations. On the 
other hand, when the Keplerian shear is included, growth 
rates higher than the shear rate 3/2I2n are expected to b e 
required to overcome the shear lllshitsu fc Sekival|2003ll . 
although numerical simulations in 3-D would be required 
to address these analytical results in detail. 

Our measured growth rates are somewhat smaller than 
in GO05, but we believe that this is due to the different 


Inp(z) = 


lnpr + ^ 

*-^s 

In Pi -h 



Fig. 14.— The logarithmic gas density as a function of height 
over the mid-plane z in the presence of falling dust. The drag force 
exerted on the gas by the falling dust mimics an extra gravity near 
the mid-plane, making the gas scale height slightly lower close to 
the mid-plane. The result is the formation of a cusp, although of 
a very moderate amplitude of about 1/10000 compared to a disk 
with no dust sedimentation (dotted line). 
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Fig. 15.— Initial growth rate of the Kelvin-Helmholtz instability 
as a function of the Richardson number Ri. There is measurable 
growth out to at least Ri = 2.0, which is way beyond the classical 
value of the critical Richardson number of Ric = 0.25. 



Coriolis force term in the present work. Changing the 
factor —1/2 to a factor —2 in equations (|3) and m 
yields very similar growth rates to GO05. The factor 
— 1/2 is a consequence of the advection of the Keplerian 
rotation velocity when fluid parcels move radially, an ef¬ 
fect that was not included in the simulations of GO05. 

7. PROPERTIES OE KELVIN-HELMHOLTZ TURBULENCE 
7.1. Diffusion Coefficient 
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The effect of turbulence on the vertical distribution 
of dust grains can be quantified as a diffu sion process 
with the turbulent diffusion coefficient Dt llCuzzLsLal] 
Il993t iDubrulle et alJlT9^ iSchraoler fc Hennind 1200 JT 
Assuming that Dt is a constant, i.e. independent of the 
height over the mid-plane, the equilibrium between ver¬ 
tical settling of dust with velocity Wz = —TfilgZ and 
turbulent diffusion implies a vertical distribution of the 
dust-t o-gas ratio efz) that is Gaussian around the mid¬ 
plane llDubrulle et al .1 ItI 9 511 . 

e = ei expl-z^/(2II^)] , (36) 

with the dust-to-gas ratio scale height given by the ex¬ 
pression iJg = Dt/(Tff2g). Writing now Dt = StH^Do, 
we get 

Doti . (37) 

Using equation ea, one can translate the scale-height of 
the dust-to-gas ratio into a turbulent diffusion coeffi¬ 
cient df Such an approach has been used to calculate the 
turbulen t diffusion coefficient of magnetorotational tur¬ 
bulence ll-Toha, risen fc Kla,hTl2nnfi|l . An obvious difference 
between Kelvin-Helmholtz turbulence and magnetorota¬ 
tional turbulence is that dust itself plays the active role 
for the first, whereas for the latter the presence of dust 
does not change the turbulence in any way, because the 
local dust-to-gas ratio is assumed to be low. Thus, for 
Kelvin-Helmholtz turbulence we expect the diffusion co¬ 
efficient to depend on the friction time 

The calculated turbulent diffusion coefficients for all 
the simulations are shown in Table|21 For the dust-to-gas 
ratio scale height we have, for simplicity, used the root- 
mean-square of the ^-coordinates of all the particles. The 
measured coefficients are extremely low, on the order of 
10“®. If we assume that the turbulent viscosity at is sim¬ 
ilar to the turbulent diffusion coefficient then one sees 
that Kelvin-Helmholtz turbulence is much weaker than 
magnetorotational turbulence where a-values from 10 “^ 
to unity are found. There is a good agreement between 
the diffusion coefficients of run B and run Br512 (which 
has twice the grid and particle resolution). This shows 
that the solution has converged and that 256 x 128 is in¬ 
deed a sufficient resolution to say something meaningful 
about the Kelvin-Helmholtz turbulence. For the simu¬ 
lations with an increased global dust-to-gas ratio (runs 
Be2, Be5 and BelO), the dust scale height falls with in¬ 
creasing dust-to-gas ratio. This is to be expected from 
equation (HOI, because of the cusp of high dust density 
that forms around the mid-plane when ei 3> 1. The 
diffusion coefficient for = 1.0 is around 50% lower 
than for fl^Tf = 0.1. Here the strong vertical settling of 
the dust has decreased the width of the dust layer sig¬ 
nificantly, and thus the diffusion coefficient is also much 
lower. 

Turbulent transport coefficients such as at and have 
an inherent dependence on the width of the turbulent 
region. Thus the “strength” of the turbulence is better 
illustrated by the actual turbulent velocity fluctuations. 
In Fig. ^1 we plot the root-mean-square of the vertical 
gas velocity as a function of height over the mid-plane. 
In the mid-plane, the value is quite independent of the 
friction time, whereas the width of the turbulent region is 
much smaller for Dgrt = 1. Thus the turbulence in itself 


TABLE 2 

Diffusion Coefficients 


Run 

17orf 


<5t/10-6 

A 

0.02 

0.0121 ± 0.0010 

3.0 ±0.5 

B 

0.10 

0.0094 ± 0.0008 

8.9 ±1.5 

C 

1.00 

0.0025 ± 0.0004 

6.5 ±2.4 

Be2 

0.10 

0.0081 ± 0.0006 

6.5 ±1.0 

Be5 

0.10 

0.0047 ± 0.0003 

2.2 ±0.3 

BelO 

0.10 

0.0031 ± 0.0002 

1.0 ±0.1 

Br512 

0.10 

0.0086 ± 0.0005 

7.5 ±0.8 


Note. — Col. (1): Name of run. Col. (2): 
Stokes number. Col. (3): Dust scale height. 
Col. (4): Diffusion coefficient derived from 
equation G3- 


is not weaker, only the turbulent region is smaller, and 
that means that the transport coefficients are accordingly 
small. 


7.2. Comparison With Analytical Work 

It is evident from Table El that the diffusion coefficient 
depends on the friction time. In the limit of small Stokes 
numbers, the constant Richardson number density dis¬ 
tribution formulated by S98 predicts that the vertical 
dust density distribution should not depend on friction 
time, and thus, according to equation that the dif¬ 
fusion coefficient should be proportional to the friction 
time. The ratio of the diffusion coefficient of run B to 
that of run A is 8 .9/3.0 ~ 3, and not 0.1/0.02 = 5 as 
would give rise to the same scale height. The factor two 
difference can be (trivially) attributed to the slight dif¬ 
ference in scale heights for the two runs. The squared 
value is different by almost a factor of two, an indica¬ 
tion that vertical settling is not completely negligible for 
= 0 . 1 . 

The strength of the Kelvin-Helm holtz turbulence ha s 
also been estimated analytically bv iCuzzi et a D (dH). 
They find that the turbulent viscosity vt should be ap¬ 
proximately (their equation [ 21 ]) 


(vvk)'^ 


(38) 


where vk is the Keplerian velocity, r] is the pressure gra¬ 
dient parameter defined in equation and Re* is the 
critical Reynolds number at which the flow becomes un¬ 
stable. This value can be approximated by the Rossby 
number Ro, the ratio bet ween the advection and Co riolis 
force terms of the flow. iDobrovolskis et ^ ll1 999f) esti¬ 
mate a value of Ro « 20 ... 30. Using the approximation 
77 « Cg/u^, equation ll^ can be written as 


vt — Dt — 2 G ^0 

Ro 


(39) 


which appears in its dimensionless form simply as 


5t = 



(40) 


With /3 = —0.1, equation © gives 77 = 0.003. Using 
Ro = 20... 30, equation (EHll gives (5t ~ 3... 8 x 10“®, 
quite comparable to the values in Table El On the 
other hand, equation does not produce a dust den¬ 
sity distribution with a constant Richardson number. 
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Fig. 16. — Root-mean-square of the vertical gas velocity as a 
function of height over the mid-plane. The value is quite indepen¬ 
dent of Stokes number, but the width of the turbulent region is 
very small for = 1 because of the strong vertical settling of 

m-sized boulders. The boundary conditions set the vertical speed 
to zero at the boundaries. 


and would thus greatly overestimate the diffusion co- 
efficient for even smaller grain s. This is also noted by 
ICuzzi fc Weidenschillind (l200fiD who constrain the valid¬ 
ity of equation P|l to iloTf > 0.01. Thus our measured 
diffu sion coefficients are actually in good agreement with 
both ICuzzi et alJ lIlOO.lIl and with S98. 

8. CONCLUSIONS 

The onset of the Kelvin-Helmholtz instability in pro¬ 
toplanetary disks has been known for decades to be the 
main obstacle for the formation of planetesimals via a 
gravitational collapse of the particle subdisk. Thus the 
study of the Kelvin-Helmholtz instability is one of the 
most intriguing problems of planetesimal formation. It 
is also a challenging problem to solve, both analytically 
and numerically, because of the coevolution of the two 
components gas and dust. Whereas turbulence normally 
arises from the gas flow alone, in Kelvin-Helmholtz tur¬ 
bulence the dust grains take the active part as the source 
of turbulence by piling up around the mid-plane and thus 
turning the energetically favored vertical rotation profile 
into an unstable shear. Planetesimal formation would be 
deceptively simple could the solids only sediment unhin¬ 
dered, but nature’s dislike of thin shear flows precludes 
this by making the mid-plane turbulent. 

In the current work we have shown numerically that 
when the dust particles are free to move relative to the 
gas, the Kelvin-Helmholtz turbulence acquires an equi¬ 
librium state where the vertical settling of the solids is 


balanced by the turbulent diffusion away from the mid¬ 
plane. For cm-sized pebbles and dm-sized rocks, we find 
that the dust component forms a layer that has a con¬ 
stant Richardson number . We t hus confirm the analyti¬ 
cal predictions bv ISekiv^ II1998I1 for the first time in nu¬ 
merical simulations. 

In the saturated turbulence we find the formation of 
highly overdense regions of solids, not in the mid-plane, 
but embedded in the turbulent flow. The clumping 
is very related to the str eaming instability found by 
lYoudin fc Goodni^ ll200,^ . Dust clumps with a den¬ 
sity that is equal to or higher than the gas density orbit 
at the Keplerian velocity, so the clumps overtake sub- 
Keplerian regions of lower dust density. Thus the dense 
clumps continue to grow in size and in mass. The fi¬ 
nal size of a dust clump is given by a balance between 
this feeding and the loss of material in a rarefaction tail 
that is formed behind the clump along the sub-Keplerian 
stream. The gravitational fragmentation of the single 
clumps into planetesimals is more likely than the whole 
dust layer fragmenting, because the local dust density 
in the clumps can be more than an order of magnitude 
higher than the azimuthally averaged mid-plane density. 
This process is very much related to the gravoturbulent 
formation of pla netesimals in turbulent magnetohy dro- 
dynamical flows ll Johansen. Klahr. fc Henningll200^ . 

A full understanding of the role of Kelvin-Helmholtz 
turbulence in protoplanetary disks must eventually rely 
on simulations that include the effect of the Keplerian 
shear, so future simulations have to be extended into 
three dimensions. One can to first order expect that 
growth rates of the KHI larger than the shear rate l7o are 
required for a mod e to grow in amplitude faster than it is 
being sheared out lllshitsu fc Sekival200.‘lll . but so far it is 
an open question in how far the radial shear changes the 
appearance of the self-sustained state of Kelvin-Helmholz 
turbulence. Including furthermore the self-gravity be¬ 
tween the dust particles, it will become feasible to study 
the formation of planetesimals in one self-consistent com¬ 
puter simulation and possibly to answer one of the out¬ 
standing questions in the planet formation process. 


Computer simulations were performed at the Danish 
Center for Scientific Computing in Odense and at the 
RIO and PIA clusters at the Rechenzentrum Garching. 
We would like to thank Andrew Youdin for his critical 
reading of the original manuscript. We are also grate¬ 
ful to Gilberto Gomez for helping to find the reason for 
the difference in the KHI growth rates between our sim¬ 
ulations and his own. Our work is supported in part 
by the European Community’s Human Potential Pro¬ 
gramme under contract HPRN-CT-2002-00308, PLAN¬ 
ETS. 
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